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An Improved Model for the Dynamical Evolution of Dark 
Matter Subhaloes 
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ABSTRACT 

Using an analytical model, we study the evolution of subhalo, including its mass, angu- 
lar momentum and merging time-scale. This model considers the dominant processes 
governing subhalo evolution, such as dynamical friction, tidal stripping and tidal heat- 
ing. We find that in order to best match the evolution of angular momentum measured 
from N-body simulation, mass stripping by tidal force should become inefficient after 
subhalo has experienced a few passages of pericenter. It is also found that the often 
used Coulomb logarithm In M/m has to be revised to best fit the merging time-scales 
from simulation. Combining the analytical model with the Extended Press-Schechter 
(EPS) based merger trees, we study the subhalo mass function, and their spatial dis- 
tribution in a Milky- Way (MW) type halo. By tuning the tidal stripping efficiency, we 
can gain a better match to the subhalo mass function from simulation. The predicted 
distribution of subhaloes is found to agree with the distribution of MW satellites, but 
is more concentrated than the simulation results. The radial distribution of subhaloes 
depends weakly on subhaloes mass at both present day and the time of accretion, but 
strongly on the accretion time. Using the improved model, we measure the second 
moment of the subhalo occupation distribution, and it agrees well with the results of 
Kravtsov et al. (2004a) and Zheng ct al. (2005). 
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1 INTRODUCTION 

In the popular cold dark matter model, structure (dark mat- 
ter halo) formation is processed in a hierarchical manner 
that small haloes form first, and they subsequently merge 
to form bigger haloes. The relics of merging haloes are seen 
as the normal galaxies in clusters, or dwarf satellites in the 
Milky- Way. In the context of galaxy formation, halo mergers 
play an important role, as they can significantly affect the 
star formation rate and morphology of galaxies. It is now 
widely accepted that elliptical galaxies are formed by major 
mergers (e.g., Toomre & Toomre 1972), and disk galaxies 
may experience preferentially minor mergers, or earlier ma- 
jor mergers if any. Thus one important aspect about galaxy 
formation in the cold dark matter (CDM) scenario is to un- 
derstand how and when the mergers (halo merger) happen, 
how the mass and density profile of accreted haloes evolve. 
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and what are their final fates: merge with central galaxies 
or get disrupted before sinking into the halo center. 

The only appropriate way to study the properties of 
accreted haloes (subhaloes) is the fully dynamic-traced sim- 
ulation. Earlier simulations (e.g., Katz & White 1993) suf- 
fer significant resolution effects, and they produce the over- 
merging pictures (e.g., Klypin et al. 1999; Moore et al. 1999). 
Simulations with higher resolution (Springel et al. 2001; Die- 
mand et al. 2004; Gao et al. 2004b; Kang et al. 2005), es- 
pecially the recent ones from two groups (Via Lactea: Die- 
mand, Kuhlen & Madau 2007a; Aquarius: Springel et al. 
2008) are shown to be capable of resolving subhalo down 
to very low mass and these simulations converged on the 
statistical distributions of subhaloes. For example, the sub- 
halo mass function (SHMF) is found to be well described by 
a single power law in both low and high-mass host haloes. 
Normalized by the host halo mass, the SHMF is universal 
with a slight dependence on formation time of the host halo 
(Gao et al. 2004b; Kang et al. 2005; van den Bosch et al. 
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2005). Those high-resolution N-body simulations also agree 
on the radial distribution of subhaloes and it is found to be 
shallower than the dark matter particles, and can be well fit- 
ted by an Einasto profile (e.g., Diemand et al. 2004; Springel 
et al. 2008). Other properties of subhaloes are also discussed 
but with diverse conclusions, such as the velocity bias of 
subhaloes, both a positive velocity bias (e.g., Diemand et 
al. 2004) and negative bias (Springel et al. 2001) are re- 
ported. The density profile of subhalo is rapidly truncated, 
with a higher concentration, but disagreements are hold for 
its inner density slope (Hayashi et al. 2003; Kazantzidis et 
al. 2004a; Diemand, Kuhlen & Madau 2007b; Springel et al. 
2008). 

In addition to the above statistical distributions of sub- 
haloes, many studies also focus on their dynamical evolu- 
tion. It is widely agreed that subhaloes will sink towards 
the host halo center by dynamical friction and gradually 
lose their mass due to tidal stripping. There are disagree- 
ment on the inner structure of subhaloes, especially at their 
late stages of evolution. Hayashi et al. (2003) found that 
subhaloes will redistribute their inner mass by tidal heat- 
ing, and become disrupted soon once the tidal radius are 
smaller than their characteristic radius. Kazantzidis et al. 
(2004a) argued that the inner part of subhalo is resistant 
to tidal shock and subhalo can orbit in the host halo for 
a longer time. They pointed out that the rapid disruption 
of subhaloes is unsurprised if the initial conditions in sim- 
ulations are not constructed in equilibrium. Kazantzidis et 
al. (2004b) further showed that numerical effects can also 
lead to rapid loss of mass from subhalo. Diemand et al. 
(2007b) found that subhaloes in their simulation can sur- 
vive for longer time even after they have passed the very 
central part of host halo where the tidal force is very strong. 
The fate of subhalo is even more complicated by the pres- 
ence of baryon. A few simulations (Gnedin et al. 2004; Nagai 
& Kravtsov 2005; Maccio et al. 2006; Weinberg et al. 2008; 
Dolag et al. 2009) have shown that compared to pure dark 
matter simulations, smooth particles hydrodynamic (SPH) 
simulations with baryon will leave more subhaloes in the 
host center as the condensation of baryon cores makes sub- 
haloes more resistant to tidal disruption, and produce a ra- 
dial distribution of subhaloes similar to that observed for 
the Milky Way satellites. 

Among the studies of subhalo dynamical evolution, one 
important issue is how long it takes for a subhalo to sink 
into the center of its host halo. This is very important for the 
model of galaxy formation as it determines when the mergers 
of galaxies actually happen. This time-scale is often called 
as the dynamical friction time scales (Tdf ). Tdt was firstly 
derived by Binney & Tremaine (1987, hereafter BT87) and 
Lacey & Cole (1993) based on the Chandrasekhar (1943) 
description. Early simulation (Navarro et al. 1995) found 
that the BT87 formula matches well with the simulation re- 
sults. Recently Jiang et al. (2008) and Boylan-Kolchin et 
al. (2008, hereafter BK08) both find that the BT87 formula 
under-estimates the merger time-scales, and they point out 
that BT87 neglected the mass loss of subhalo during it evo- 
lution. But even the results of Jiang et al. (2008) and BK08 
differ by a factor of two, and this diversity is from various 
effects. Jiang et al. (2008) use cosmological simulation with 
star formation and feedback, while BK08 use controlled two- 



halo merging simulation with pure dark matter. Also they 
adopt different definitions for galaxy/subhalo mergers. 

Although numerical simulation is the only proper 
method to study the dynamical evolution of subhalo, use- 
ful insight into physical processes governing subhalo evolu- 
tion can be gained from analytical model. Base on the pi- 
oneer work of Taylor & Babul (2001), the analytical model 
was well developed in the past years (Benson et al. 2002; 
van den Bosch et al. 2005; Taylor & Babul 2004; 2005a; 
2005b; Zentner & Bullock 2003; Zentner et al. 2005). The 
model includes the main physical processes governing sub- 
halo evolution: gravitational force, dynamical friction, tidal 
stripping, tidal heating and disruption. Coupled with merger 
trees from the extended Press-Schechter (EPS) theory, the 
analytical model is capable of producing realistic catalogue 
of subhaloes in given host halo, which can be directly com- 
pared to N-body simulation results. Up to now most of these 
analytical works unfortunately neglect the study of subhalo 
merging time-scales. Although Taffoni et al. (2003) derived 
an fitting formula for Tdf , their results are not tested against 
simulation results. BK08 recently found that their results 
are still quantitatively inconsistent with the prediction of 
Taffoni et al. (2003). 

In this paper, using an analytical model similar to Zent- 
ner et al. (2005; hereafter Z05), we study the dynamical evo- 
lution of subhalo. We investigate the effects of tidal strip- 
ping. Coulomb logarithm on the angular momentunij evolu- 
tion, merging time-scale of subhalo. By comparing our model 
predictions to the simulation results of BK08, we proposed 
a modified Coulomb logarithm which can well reproduce the 
evolution of angular momentum and merging time-scales for 
subhalo with different mass ratio and orbit eccentricity. We 
then combine the analytical model with the Monte Carlo 
merger tree to produce subhalo catalogue in a Milky Way 
(MW) type halo, and compare the model predictions to both 
N-body simulation results and observed distribution of satel- 
lites in the Milky Way. In Section [2] we present the main 
ingredients of the analytical model and show the model pre- 
dictions in Section O In Section U we further examine our 
model with the halo occupation distribution of subhaloes. 
In section [5] we discuss the radial distribution of subhaloes, 
and we briefly conclude our model in Section [S] 

Throughout this paper we adopt a fiat ACDM cosmol- 
ogy with the following cosmological parameters: Q.m = 0.25, 
Oa = 0.75, h = //o/(100 kms"^ Mpc"^) = 0.73, fib = 0.04, 
Ws = 0.951 and erg = 0.9. 



2 THE MODEL 

In this section we present the model for the evolution of 
the population of dark matter subhaloes. The first part de- 
scribes the merger history (i.e. mass assembly) of the host 
halo, which can be obtained using either the EPS theory or 
TV-body simulations. The second part describes the dynam- 
ical evolution of subhaloes, and includes orbit integration 
in the presence of dynamical friction combined with tidal 
stripping and heating. Our model is similar to that of Z05, 

^ Throughout this paper, when referring to the angular momen- 
tum, we mean it is the angular momentum per unit mass or the 
specific angular momentum. 
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Figure 1. The 10 random mass accretion histories (MAHs, thin 
lines) for the mass of the main haloes as function of time, for a 
MW sized halo with M(z = 0) = 1.77 X IQ^'^h-^MQ. The thick 
line with error bars show the average MAH and 1 — a scatters 
from 100 such reahzations. 



but there are also a few differences. We employ the well- 
calibrated code of Parkinson et al. (2008) to construct the 
EPS merger trees. This code is a significant improvement 
over the Somerville & Kolatt (1999) implementation used 
by Z05. In addition, we calibrate our model using detailed 
numerical simulations on the evolution of the orbital angular 
momentum of subhaloes. As we will demonstrate, this kind 
of calibration is far more constraining than using the mass 
function or velocity function of subhaloes. Finally, we (i) use 
a more detailed, empirical treatment of tidal heating, based 
on the results of high resolution numerical simulations, (ii) 
investigate the impact of changes in the Coulomb logarithm 
used in the analytical treatment of dynamical friction, and 
(iii) consider a different treatment for the tidal mass loss of 
subhaloes. 

In what follows, we use m and M to denote the in- 
stantaneous masses of subhalo and host halo, respectively. 
Unless stated otherwise, we consider it understood that both 
m and M are functions of time. Then we use the sym- 
bol /i to refer to the mass ratio between subhalo and host 
halo, i.e., /x = iJ.{t) = m{t)/M{t), without writing the time- 
dependence explicitly. We use /^i, /if to refer to the initial 
mass ratio at the time of accretion (tacc or Zacc) and the final 
mass ratio at present day [z — 0), i.e., Hi = m (face) /A/ (face), 
fit = m{z — 0)/M{z — 0), respectively. 



2.1 Merger Trees 

The backbone for modeling the population of dark matter 
subhaloes is the merger history of the host halo, which de- 
scribes when each subhalo is accreted, and what its mass 
is at accretion. Halo merger histories can be obtained using 
A''-body simulations, or in a semi-analytical fashion from 
the EPS formalism (Bond et al. 1991; Lacey & Cole 1993). 
Here we adopt the latter approach by employing the open- 
source code of Parkinson et al. (2008) to generate the as- 



sembly histories of MW sized haloes with a. z = mass of 
M{z = 0) = 1.77 X 10^^ A/0 that is close to the z = mass of 
the host halo in "Via Lactea" and "Aquarius" simulations. 
We construct 100 independent merger-tree realizations, each 
with a resolution of 10* A/©. In Fig. [l] we show an example 
of the halo mass accretion history (MAH) , which is defined 
as the mass of the most massive progenitors at each redshift 
from the merger tree. The thin lines are 10 random MAHs of 
the MW type halo in our adopted cosmology. The thick line 
with error bars is the average MAH and its 1 ~ a scatters 
(standard deviation) from 100 such realizations. 

The main branch of the merger tree is defined as the 
trajectory of the most massive progenitors starting from the 
z — halo. In our study we consider only those haloes that 
are directly accreted onto the main branch, not account- 
ing for any of their subhaloes (which would give rise to 
sub-subhaloes). As shown by Yang, Mo & van den Bosch 
(2009), sub-subhaloes (and higher-order substructures) only 
contribute a small fraction to the total substructure mass 
function (see also Giocoli et al. 2010). 

For each halo in the merger trees, its virial radius, rvir, 
is defined as the radius within which the mean mass density 
is Ac (2) times the critical density of the universe at redshift 
2, where 



Ac(2) = 187r^ + 822: - 39a;^ 



(1) 



with X — Sim (2) — 1 (Bryan & Norman 1998). We assume 
that the host haloes are spheres with a density distribution 
given by the NFW profile (Navarro et al. 1997). The corre- 
sponding concentration parameter, c, is set using the median 
relation between c and halo mass M of Neto et al. (2007), 
which is given by 



c{M,z) 



4.67 

l + z 



M{z) 



10i4/i-iM( 







(2) 



where the dependence on redshift is taken from Bullock et 
al. (2001). Dark matter subhaloes are assumed to have a 
similar NFW profile at their time of accretion (i.e., at the 
time they become a subhalo), but as described below (see 
Section [2.2.5(1 . this density profile is subsequently modified 
due to tidal heating. 



2.2 Dynamical Evolution of Subhaloes 

2.2.1 Orbital Parameters 

The first step for the dynamical evolution of dark matter 
subhaloes is to assign their orbital parameters at the time 
of accretion, face- We follow Z05 and draw the initial orbital 
energy and angular momentum from distributions that have 
been obtained from A^-body simulations. In particular, we 
assume that each subhalo starts its orbit at the virial radius, 
rvir, of the host halo at the time of accretion with an orbital 
energy equal to that of a circular orbit of radius ?7rvir, where 
rj is drawn randomly from a uniform distribution between 
[0.6, 1.0] (see Z05). The initial specific angular momentum is 
parameterized as jinit = eje, where jc is the specific angular 
momentum of the circular orbit mentioned above and e is 
called the orbital circularity (note that < e < 1). Several 
studies (e.g., Benson 2005; Tormen 1997; Z05; Khochfar & 
Burkert 2006; Jiang et al. 2008) have measured the distri- 
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bution of £, all reporting similar results. Here we use the 
distribution obtained by Jiang et al. (2008): 

\2.99 



/(e) = 2.77e''''' (1.55 - 



2.2.2 Dynamical Friction 



(3) 



We treat subhalo as test particle in the orbital evolution. In 
addition to the (radial) force due to the gravitational poten- 
tial of the spherical NFW host halo, subhalo experience an 
effective force due to 'dynamical friction' caused by the grav- 
itational interaction between subhalo and the background 
'field' particles that make up the host halo. Chandrasekhar 
(1943) showed that if the distribution of background par- 
ticles is infinite and homogeneous, one can obtain an ana- 
lytical expression for the dynamical friction force by consid- 
ering the cumulative effect of many uncorrelated two-body 
interactionqj between the subject mass (in our case the sub- 
halo) and the individual field particles. This is known as the 
Chandrasekhar dynamical friction force, which is given by 



(4) 



Fdf = -4^ (— y lnAp« i;orb)^^ . 

Here m and Voa are the mas^ and orbital velocity of the 
subhalo, In A is the Coulomb logarithm, and p(< «orb) is the 
density of the particles in the host halo that have a speed 
less than the velocity of the subhalo. The Coulomb loga- 
rithm is introduced to avoid divergence that arises from the 
assumption of an infinite, homogeneous sea of field particles. 
Similar to the frictional drag in fiuid mechanics, dynam- 
ical friction exerts a force always opposite to the motion. 
However, contrary to hydrodynamic friction, which always 
increases in strength when the velocity increases, the drag 
due to dynamical friction has a more complicated depen- 
dence on velocity. While _Fdf oc iiorb in the low iiorb-limit, 
similar to hydrodynamic friction, one has that _Fdf oc u"^ 
in the high Worb-hmit. In what follows we assume that the 
'field' particles that make up the host halo follow a locally 
Maxwellian velocity distribution. In that case. Equation Q 
reduces to 



fGmy 



Fdf = -47r ( ^^^ ) InAp(r) 



crf(A:) 



2X ^^. 



Vorb 
Worb 



(BT87), where X = Vovh /[V2cr{r)] with a{r) the local, one- 
dimensional velocity dispersion of the host halo at radius r, 
which can be solved using the Jeans equation (BT87; Cole 
& Lacey 1996) under the assumption that the stress tensor 
is isotropic]. 



2.2.3 Coulomb Logarithm 

As shown by White (1976), in the case of an extended sub- 
ject mass, as is the case for our subhalo, one has that 



^ By considering the interactions to be uncorrelated, one effec- 
tively ignores the self-gravity of the field particles. 
^ Note that the mass entering to the dynamical friction may not 
be the same as the bound mass of subhalo (e.g., Fellhauer &; Lin 
2007). Here, an 'effective' mass contributing to the dynamical 
friction is modeled. 
^ see Zentner & Bullock (2003) for a useful fitting function. 
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m(r)dr 
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(5) d X _ 



with m{r) the subhalo mass profile. Here fcmax is the max- 
imum impact parameter considered, which is introduced in 
order to avoid divergence. This divergence, however, arises 
because Equation Q is based on the (unrealistic) assump- 
tion of a homogeneous and infinite medium. In our case of a 
subhalo orbiting in a host halo, a logical value for fcmax rnay 
appear to be the size of the host system, which is indeed 
what is often adopted. However, it is important to realize 
that strictly speaking Equation Q is not valid for this case, 
and that there is no 'correct' value for the Coulomb loga- 
rithm. Hence, different forms for In A have been adopted in 
the literature. Some authors treat the Coulomb logarithm as 
a constant (Velazquez & White 1999; Taylor & Babul 2001; 
Jardel & Sellwood 2009) . Others claims that this yields a dy- 
namical friction time, Tdf, defined as the timescale on which 
the subject mass looses its orbital angular momentum, that 
is too short, and advocate instead that In A has to be time- 
dependent (e.g., Colpi et al. 1999; Hashimoto et al. 2003). A 
widely used form for the Coulomb logarithm is \n{M/m) or 
ln(l -I- M/m), where M and m are the instantaneous (time- 
dependent) mass of the host and subhalo, respectively (e.g. 
BK08; Jiang et al. 2008). In this paper we will consider two 
different forms for the Coulomb logarithm: InA = C and 
In A = — In p -I- C. As we will show, both forms yield equally 
satisfactory results (when compared to numerical simula- 
tions), as long as C is allowed to vary with the initial orbit 
parameters, r/ and e, and the initial (i.e., at the time of infall) 
mass ratio ^i. 

2.2.4 Orbit Integration 

We integrate the orbits [x(r-, 0)] of subhalo by treating it as 
test particle. The equation of motion for a subhalo of mass 
m is given by: 

d^ 
df2 



GM{< r) r Fdf 



(8) 



with _A/(< r) the mass of the host halo inside of radius r, 
and Fdf the dynamical friction force given by Equation ([5]). 
The equation of motion is solved using a fifth-order Cash- 
Karp Runga-Kutta method. During time-steps in which the 
mass of the host halo increases (due to the accretion of a new 
subhalo), we recompute the mass distribution and potential 
of the host halo, always under the assumption that the host 
halo has a NFW shape with the c{M) relation as given by 
Equation ([2]). We assume that the orbital angular momen- 
tum of a subhalo is conserved when the mass of the host 
halo increases; the only mechanism by which the subhalo 
is assumed to lose orbital angular momentum is dynamical 
friction. 



2.2.5 Tidal Stripping and Heating 

When a subhalo orbits its host halo, it loses mass due to 
tidal stripping. The tidal radius, rt, is the radius in sub- 
halo where the external differential (tidal) force from the 
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host halo exceeds the binding force of the subhalo, and is 
approximated by 



Gm{< n) 



uj^ + G [2M{< r)/r^ - inp{r)] ' 



(9) 



with cj the angular speed of the subhalo and p{r) the density 
profile of the host halo (von Hoerner 1957; King 1962; Tay- 
lor & Babul 2001). The subhalo mass outside n becomes 
unbound and is ultimately stripped. It should be pointed 
out, however, that Equation (|9]) is only a crude approxima- 
tion. First of all, in the case of non-circular orbits the con- 
cept of a tidal radius is not well defined. Secondly, even in 
the case of point masses, the two-dimensional surface along 
which d^rt/di^ — (i.e., zero- velocity surface; BT87) is not 
spherical, and so cannot be characterized by a single radius. 
And finally. Equation ((9} ignores the orbital motion of par- 
ticles within the subject mass. This, among other effects, 
gives rise to scatter in uj, and effectively introduces some 
non-zero 'thickness' to the shell of particles for which the 
internal and tidal forces balance. 

Because of these uncertainties, numerical simulations 
show that many particles remain bound even though they 
lie beyond the tidal limit of Equation ® when the subhalo 
is near pericenter (e.g., Diemand et al. 2007b). Fellhauer & 
Lin (2007) also showed that the previously stripped material 
can contribute to the dynamical friction and affect mass loss 
from the subhalo. This has resulted in uncertainties regard- 
ing how best to model tidal stripping. In particular, difi'erent 
studies adopt different time scales for tidal stripping, Tstrip, 
defined by 

dm _ m{> rt) 

Ot _£ strip 

Whereas Taylor & Babul (2001) simply assumed that Tstrip 
is equal to the instantaneous orbital time Torb = 2-k jio, Z05 
and Diemand et al. (2007b) inferred stripping time-scales 
that are 3.5 and 6 times shorter, respectively. In order to 
parameterize this uncertainty we adopt 



(10) 



Te 



To, 



strip 



(11) 



which we use in combination with Equations. ((9)1 and (|10p 
to describe mass loss due to tidal stripping. Here A is the 
tidal stripping efficiency parameter, which we tune using 
detailed numerical simulations (see Section FS. 21 below V Note 
that Taylor & Babul (2001), Z05 and Diemand et al. (2007b) 
used or advocated A ~ 1.0, 3.5 and 6.0, respectively. 

Using numerical A'^-body simulations, it has been shown 
that tidal heating causes subhaloes to expand and to reduce 
their inner mass profile (e.g., Hayashi et al. 2003; Kravtsov 
et al. 2004b). Hayashi et al. (2003) introduced a modified 
NFW profile to describe the density distribution of a tidally 
heated subhalo according to 



p{r) = 



/t 



1 + (r/rte)3 



PNFw(r) . 



(12) 



Here pnfw('') is the original NFW density profile of the 
subhalo at the time of infall, rtc is the 'effective' tidal radius 
that describes the outer cutoff imposed by the tides, and /t 
describes the reduction in the central density of the subhalo. 
As shown by Hayashi et al. (2003), these are well- fit by 



and 

Ig /t = -0.007 + 0.35 Qra + 0.39 Ql, -\- 0.23 QL • (14) 

Here Qrn = lg[?Ti(t)/"T-(tacc)] is the logarithm of the remain- 
ing fraction of subhalo mass, and r^ is the scale radius of 
the NFW profile at the time of accretion. Both /t and rtc 
decrease with time while a subhalo is losing mass. 

We caution that this 'empirical' treatment of tidal heat- 
ing is subject to some debate. In particular, Kazantzidis et 
al. (2004a) have argued that the simulation used by Hayashi 
et al. (2003) was not set-up in equilibrium, and that this 
has resulted in a tidal mass loss rate that is too high. Un- 
fortunately, lacking a more reliable description of how tidal 
heating impacts the density profiles of subhaloes, we use 
Equation (I12|l despite these potential problems. 



w !i£ = 1.02 



1.38 Qn. + 0.37 Q„ 



(13) 



2,. 2. 6 Disruption and Cannibalism 

As a subhalo is being exposed to tidal stripping and heating, 
it may reach a point at which it is disrupted, i.e., at which no 
significant amount of matter remains gravitationally bound 
into a single object. Alternately, depending on the dynamical 
friction time, a subhalo may either sink all the way to the 
center of the host halo's potential well (i.e., lose all its orbital 
angular momentum) or continue to orbit as a subhalo, if 
the mass ratio m{t)/M{t) is such that dynamical friction is 
negligible. 

Whether and when subhaloes are tidally disrupted is 
still being debated. Testing this with numerical simulations 
is complicated by the fact that simulations are always sub- 
ject to numerical artifacts due to limited mass and force 
resolution. Using A'^-body simulations, Hayashi et al. (2003) 
found that subhaloes are disrupted once their tidal radius 
rt becomes smaller than ~ 2rs. Motivated by these findings, 
Taylor & Babul (2004) and Z05 included tidal disruption in 
their semi-analytical models. They considered a subhalo to 
become tidally disrupted once its mass becomes less than 
its initial mass (i.e. at accretion) within a radius /disrs, with 
/dis = 0.1 (Taylor & Babul 2004) and /dis = 1.0 (Z05), 
respectively. Using the luminosity function of Milky Way 
satellites, Maccio et al. (2009) argued that 0.1 < /dis ^ 0.5, 
while Wetzel & White (2010), using a wide variety of obser- 
vational constraints on satellite galaxies, conclude that sub- 
haloes are disrupted once their bound mass drops below ~ 2 
percent of its mass at infall, corresponding to /dis ~ 0.3 for 
a NFW profile with a concentration c = 10. Hence, typical 
values for /dis that have been adopted in the literature cover 
the entire range 0.1 <, /dis ^ 2; as nicely shown by Wetzel 
& White (2010), varying /dis by this amount has a huge im- 
pact on the radial number density distribution of surviving 
subhaloes, with smaller /dis resulting in a more concentrated 
profile. 

There have also been a number of studies that have ar- 
gued that subhalo disruption is actually extremely rare. In 
particular, a number of numerical simulations with very high 
spatial resolution have shown that subhaloes are remarkably 
resilient to disruption by tidal shocks (e.g., Kazantzidis et 
al. 2004a; Bullock & Johnston 2005; Penarrubia et al. 2008; 
Springel et al. 2008). Kanzantzidis et al. (2004a) argued 
that the initial conditions of the simulation of Hayashi et 
al. (2003) were not in equilibrium, which is likely to have 
caused a subhalo disruption rate that is too high. Using the 
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"Via Lactea" simulation, Diemand et al. (2007b) found that 
only very few subhaloes get disrupted; In their simulation 
97% of the subhaloes identified at z — 1 were still present 
at z = 0, and even subhaloes with rt < 0.2rs were found to 
survive. Finally, BK08 also found that subhaloes survive as 
bound entities up to the point of having lost all their orbital 
angular momentum (private communication). 

Motivated by these high-resolution simulations, we as- 
sume that subhaloes are never disrupted. Rather, when a 
subhalo has lost all its orbital angular momentum, we con- 
sider it 'cannibalized' by (or 'merged' with) the host halo, 
and we remove the subhalo from our samplqj. The merg- 
ing time scale, Tdf, is defined as the time interval between 
accretion and merging of a subhalo. We consider two cases 
for subhalo mass loss due to tidal stripping. In the first case 
(hereafter Model Ml), subhaloes can lose mass continuously 
until they are cannibalized by their host halo. In the second 
case (hereafter Model M2), tidal stripping is 'turned off' 
after a subhalo has experienced two pericentric passages. 
This is motivated by numerical simulations, which suggest 
that subhaloes only experience significant mass loss during 
their first two orbital periods (e.g., Diemand et al. 2007b). 
As we will show below (see Section [3. 1|) . these two different 
treatments of tidal stripping yield very different predictions 
regarding the evolution of the orbital angular momentum of 
subhaloes. 



3 COMPARISON WITH NUMERICAL 
SIMULATIONS 

We now turn to a detailed comparison of our analytical 
model with numerical simulations. After exploring how the 
orbital evolution of a subhalo depends on the tidal strip- 
ping efficiency. A, and the Coulomb logarithm. In A, we tune 
these parameters by fitting the merging time-scales and the 
evolution of orbital angular momentum of subhaloes to the 
controlled, high-resolution numerical simulations of BK08. 
These simulations follow the orbital evolution of individual 
subhaloes of different mass and with different orbital prop- 
erties in a host halo of fixed mass M{t) = M, and are ideally 
suited to tune our model parameters. 

Subsequently, we use our model to compute the mass 
function and radial number density distribution of subhaloes 
in a MW type host halo, which we compare to numerical sim- 
ulations and observational constraints from satellite galaxies 
in the MW. 



3.1 Tuning the Tidal Stripping Efficiency and 
Coulomb Logarithm 

Fig. [5] shows the evolution of the mass (upper panels) and 
halo-centric radius (lower panels; in units of the virial radius 
of the host halo) of a subhalo with /ii = 0.05, e = 0.5 and rj — 
1.0 in the model Ml. In the left-hand panels the Coulomb 
logarithm is In A = — In /x and we vary the efficiency of tidal 

^ Note that some authors consider subhaloes 'merged' or 'canni- 
balized' when its separation from the host center is smaller than 
some fiducial radius (e.g., Kravtsov et al. 2004a; Z05). However, 
we consider our definition, based on the complete loss of orbital 
angular momentum, more realistic (see also BK08). 
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Figure 3. The evolution of angular momentum with /li = 0.1, 
£ = 0.65, rj = 1.0, InA = -In^ from model Ml. The solid 
lines are the model results with different tidal stripping efficiency 
{A = 1, 3.5, 10). The dashed line shows the simulation result from 
BK08. 



stripping, expressed by the parameter A (see Equation [TT]), 
as indicated. Clearly, for larger values of A (i.e., more rapid 
mass loss), the effect of dynamical friction is reduced, and 
the orbital decay slows down drastically. With decreasing 
dynamical friction, subhalo can also travels to higher halo- 
centric radius (lower panel). 

In the right-hand panels, we keep A fixed as 3.5 and 
vary the Coulomb logarithm, as indicated. Typically, a lower 
value for InA results in dynamical friction being less effi- 
cient (cf. Equation [5]). This in turn implies a reduced mass 
loss, because the subhalo experiences fewer pericentric pas- 
sages and, on average, more orbits at larger halo-centric radii 
where the tidal forces due to the host halo are weaker. 

We now turn to a detailed comparison with the sim- 
ulation results of BK08. To that extent, we use the same 
initial conditions, such as the density profiles of subhalo 
and host halo, orbital circularity and orbital energy. The 
dashed line in Fig. [3] shows the evolution of j{t) for a sub- 
halo with Hi = 0.1, e — 0.65, and rj — 1.0 in the simulation of 
BK08. The solid lines show the predictions from our model 
Ml, for three different values oi A, as indicated. In all three 
cases we have used In A = — In /x. Clearly, the evolution of 
orbital angular momentum is a strong function of the ef- 
ficiency of tidal stripping, A. Larger values of A result in 
lower rate of angular momentum loss (i.e., a longer merging 
time scale, Tdt), simply because a less massive subhalo expe- 
riences weaker dynamical friction. Whereas our model with 
A — 3.5 matches j{t) in the BK08 simulation reasonably 
well for the first ~ 5 Gyr, the predicted evolution in orbital 
angular momentum at later stages is too weak. We have ex- 
perimented with different values of A but where unable to 
obtain a satisfactory match to the j{t) in the simulation of 
BK08. 

In Fig. 21 we investigate the impact of changing the 
Coulomb logarithm. In the left-hand panels we keep InA 
fixed at some constant values (as indicated), while in the 
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Figure 2. The evolution of mass and lialo-centric radius for a sublialo witii ^i = 0.05, e = 0.5, and r) = 1.0 in model Ml. In the left-hand 
panels, In A = — In^ and we vary the efficiency A of tidal stripping, as indicated. In the right-hand panels, A = 3.5 and we vary the 
value of the Coulomb logarithm, as indicated. 
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Figure 4. As Fig. |3] but with A is fixed as 3.5 (constrained 
from subhalo mass function, shown in Fig. [Sj . The figure shows 
the effects of varying Coulomb logarithm. Note the predicted slow 
evolution at later stages. 



right-hand panels we adopt \nA = —\nfi + C for three dif- 
ferent values of C (as indicated). In all case A is fixed as 3.5, 
which, as we will see in Section [3.21 yields a subhalo mass 
function that is in best agreement with numerical simula- 
tions. Although different Coulomb logarithms have a signifi- 
cant impact on the evolution of the orbital angular momen- 
tum, we were unable to find a form for In A for which we 
could satisfactorily match the simulation results of BK08, 
even if we kept A a free parameter as well. The problem 



is that the model typically predicts a decline in the angu- 
lar momentum loss rate, 7?. = —dj/dt, while the simulation 
results have 7?. ~ constant during the entire evolution, up 
to the point of being cannibalized by the host (i.e., when 
j = 0). The only exception is when the tidal stripping effi- 
ciency A is very low, in which case the merging time scale, 
Tdf , is much too short. The culprit for this discrepancy is the 
continued mass loss due to tidal stripping. This has moti- 
vated us to consider a modified model (model M2) in which 
tidal stripping is inefficient after two pericentric passages. 
As already mentioned above, this is not an entirely ad-hoc 
modification, as it has support from the ultra-high resolu- 
tion "Via Lactea" simulation (Diemand et al. 2007b). 

Fig. [5] shows the evolution of the orbital angular mo- 
mentum of two different subhaloes from the model M2; in 
the upper panels /ii — 0.1, e = 0.65, and i] — 1.0, while the 
lower panels correspond to a subhalo with the same rj but 
with fii = 0.05 and e = 0.46. The dashed lines indicate the 
results from the simulations of BK08, while the solid lines 
correspond to our model M2 with A = 3.5, and with the 
Coulomb logarithm tuned to best match the BK08 results. 
In the left- and right-hand panels we considered InA = C 
and In A = — In fi + C, respectively, where C is a constant. In 
all four cases, the fit to the j{t) of BK08 is fairly satisfactory. 
Clearly, model M2 gives much better fit to the simulation 
results, and we adopt this model throughout the remainder 
of this paper (unless specifically stated otherwise). 
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Figure 5. The evolution of angular momentum in model M2, in 
which tidal stripping is stopped after subhalo has gone through 
two pericenter passages. The upper panels are evolutions with 
/i; = 0.1, £ = 0.65 and tj = 1.0. Lower panels are for subhalo 
with ^i = 0.05, e = 0.46, rj = 1.0. The two forms of Coulomb 
logarithm that In A = C and In A = — In /.t + C are used in left 
and right panels, respectively. The solid and dashed lines indicate 
the model and simulation {BK08) predictions. Better agreement 
between them can be obtained with appropriate In A. 

Note, though, that the best-fit value of C (indicated in 
each panel) is different in each case. Both the model with 
in A = C and that with In A = — In ^ + C yield equally 
satisfactory results. In what follows we will only consider 
the latter, since we believe it to be the more physical one. 
What remains to be done, however, is to characterize how C 
depends on the mass and orbital properties of the subhalo. 
Using a suite of numerical simulations, BK08 derived a fit- 
ting formula for the merging time scale, Tdf, as a function of 
the mass, m, the orbital circularity e, and the orbital energy, 
T], of the subhalo at accretion. We use this fitting formula to 
constrain C = C[/ii,?7,e]. After some experimenting, we fi- 
nally adopted the following form for the Coulomb logarithm: 



In A = — In/i -f CYi^^jf^ exp[c4e] -I- C5 . 



(15) 



Fitting the merging time scales in our semi-analytical model 
to those listed in Table 1 of BK08, we obtain; c\ — 1.04, 
C2 = -0.64, C3 = 0.72, C4 = -3.02, and cs = -0.75. 

Fig. [6] shows the merging time-scales as function of the 
initial mass ratio, Hi, for three different values of the orbital 
circularity, as indicated in each panel. In all cases the initial 
orbital energy has rj = 1.0, and we have adopted A = 3.5. 
The solid lines correspond to the predictions of model M2 
using the above Coulomb logarithm of Equation p5|l , while 
the dashed lines indicate the fitting formula of BK08. Al- 
though not perfect, our model is in fair agreement with the 
simulation results of BK08. This is also evident from Fig. [T] 
where we show the evolution of orbital angular momentum 
for sbc different combinations of pi and e (as indicated). In 
each panel the dashed curve corresponds to the simulation 
results of BK08, while the solid line is our M2 model predic- 
tion based on the Coulomb logarithm of Equation (|15p . The 
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Figure 7. The evolution of angular momentum for subhalo. The 
dashed lines are from simulation by BK08, and solid lines are 
model predictions. Here the Coulomb logarithm is from Equa- 
tion I I15II . and C is labeled in each panel. It shows that our revised 
Coulomb logarithm can well describe the evolution of angular mo- 
mentum for subhaloes with different mass ratio, orbit circularity. 
T] = 1.0 in all cases. 



corresponding value of C is indicated in each panel. Overall, 
our model yields j{t) that are in satisfactory agreement with 
the simulation results of BK08. 



3.2 The Distribution of Subhalo Population 



In Section 12.21 we have introduced in detail the model for 
the evolution of subhalo, including its mass, radial position 
and merging time-scales. In Section TS. II we tune the model 
parameters to fit the dynamical evolution of subhalo pre- 
dicted by simulations. Couple with the merger trees, the 
model is ready to produce the subhalo catalogue in the host 
halo. As described in Section 12.11 our model employs 100 
realizations of merger trees of the MW type halo, and each 
realization specifies a random assembly history of dark mat- 
ter haloes (Fig.[l|. We follow the dynamical evolution of the 
accreted subhaloes [with masses ?n(2acc) > 10* M©] by the 
main branch of merger tree, and investigate the distribution 
of subhalo population at z = 0, including the subhalo mass 
function (SHMF) and their radial distribution. We also com- 
pare the model prediction with the simulation results and 
the observed distribution of the MW satellites. 

Fig|S]show the SHMF and the radial number distribu- 
tion in the upper and lower panels, respectively. For panels 
in each column, the same set of model parameters is used 
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Figure 6. Subhalo merging time-scales. As in BK08, subhalo merges with central halo once it loses all of its angular momentum j. 
The three panels are for different orbital circularity e = 0.3,0.5,0.7. The solid lines are from model M2 with r) = 1.0, A = 3.5 and 
In A = — In/x + C given by Equation I I15II . The dashed line shows the fitting formula from simulations of BK08. 
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Figure 8. Subhalo mass function (SHMF, upper panels) and radial number distribution (lower panels) in a Milky- Way type halo. Panels 
a and d basically show the predictions from model M2 (solid) with A = 3.5 and Coulomb logarithm of Equation l|15|l . In panel a, the 
unevolved and evolved SHMFs from simulation are shown as dashed-dotted and long dashed lines, respectively. In panel d, we also plot 
the radial number distribution of simulated subhaloes (hatch area, upper limit: Via Lactea; lower limit; Aquarius) and observational 
MW satellites (squares). In panels b and e, we compare the results of simulation to the model Ml with In A = — In /i and A = 1, 3.5, 10, 
while in panels c and f, the used parameters are In A = 3,4, 5 and A = 3.5. The model predictions with different parameters are plotted 
in lines with varying line style as indicated. 



10 Gan et al. 



and indicated in the lower panel. The SHMF from N-body 
simulations is well described by a power law, with index be- 
tween -0.8 and -1.0 (Springel et al. 2001; Gao et al. 2004b; 
Kang et al. 2005; Diemand et al. 2004; Giocoh et al. 2008; 
2010). In panel a, we show the simulated one from Giocoli et 
al. (2008) as the long-dashed line and the model prediction 
is shown as the solid line. It can be seen that our fiducial 
model (Model M2 with A — 3.5 and Coulomb logarithm 
from Equation <\15^ ) produces a fair match to the simula- 
tion result. 

As the SHMF is for an evolved population of accreted 
subhaloes, it is important to check if the unevolved SHMF, 
which is the mass function of subhaloes at their accretion 
times, is reproduced by the EPS based merger tree employed 
in our model. The unevolved SHMFs from simulation and 
the EPS model are shown as dashed-dotted and dotted lines, 
respectively. Their good agreement indicates that the model 
for the dynamical evolution of subhalo is not biased by the 
formation history of the host halo. 

We further check if our model predictions are affected 
by the assumptions for the dynamical processes of subhalo. 
Panel b and c show the predictions from our Model Ml, 
with the dependence on A (panel b) and Coulomb logarithm 
(panel c). It is found that the SHMF depends strongly on 
tidal stripping efficiency A, but weakly on Coulomb loga- 
rithm. This can be understood from that, as shown by van 
den Bosch et al. (2005), the subhaloes population at present 
day is dominated by the recent (the last ~ 1 — 2 Gyr) accre- 
tion history of the host halo. It is already shown in Fig. [2] 
that subhalo mass depends strongly on A at the first few 
Gyrs, but with a weak dependence on Coulomb logarithm. 
Thus the results indicate that SHMF can not be used to con- 
strain the mass evolution of subhalo after a few Gyrs, while 
the dynamical evolution j can set strong constraints on the 
late stage evolution of subhalo, as shown in Section [3. II 

The radial number distribution of subhaloes is shown in 
the lower panels of Fig|8] In panel d, the hatched area shows 
the spanned distribution from simulations (upper limit: Via 
Lactea; lower limit: Aquarius). The observed distribution of 
the MW satellites is shown as the empty squares (data are 
from Mateo 1998; Kroupa et al. 2005; Metz et al. 2007; Metz 
et al. 2009; Martin et al. 2008). A clear discrepancy is that 
the distribution of the MW satellites is more concentrated 
than the subhaloes from N-body simulations, which has al- 
ready been noted before (e.g. Taylor et al. 2005b). Such a 
discrepancy could be due to the incompleteness of observa- 
tions ( Willman et al. 2004) , or the observed satellites present 
a biased population of subhaloes from simulations (Kravtsov 
et al. 2004b; Madau et al. 2008). We leave more discussion 
to Section [5l 

The fiducial model prediction is shown as the solid line 
in panel d. Compared to the simulation result, the model 
predicts a more centrally concentrated distribution of sub- 
haloes. A Similar discrepancy was also noted by Taylor & 
Babul (2005b) although their model prediction is slightly 
lower than ours. However, Z05 found that their model pre- 
dicts a well match to simulation result, and they argued 
that the discrepancy noted by Taylor & Babul (2005b) is 
not from numerical effects of simulation but the model as- 
sumptions for subhalo merging and disruption. There still 
lacks detailed studies on this issue. Here we firstly explore 
if the predicted distribution of subhaloes is affected by the 



model assumption. The predictions from our Model Ml are 
shown in panel e and panel f, with dependence on A and 
Coulomb logarithm, respectively. Surprisingly, we find that 
the predicted distribution is similar to that obtained from 
our Model M2, and it also has no dependence on the model 
parameter A and In A. 

In principle, the final spatial distribution of subhaloes 
is mainly determined by (1) their initial positions at accre- 
tion, (2) dynamical processes governing subhalo evolution, 
and (3) criteria on where and when subhalo disappears. The 
results in panel e and f suggest that (2) has no significant 
effects on the radial distribution of subhaloes. Since the low- 
mass subhaloes dominate the subhaloes population, varying 
the strength of the dynamical friction and tidal stripping 
will not change the spatial distribution of subhaloes much. 
In addition, Kang (2008) has shown that the formation his- 
tory of the host halo from the EPS theory is very similar 
to that of simulations. As the mass and radius of halo are 
close related by Equation ([l|, the initial positions of sub- 
haloes at accretion (the virial radius of host halo) from the 
EPS merger tree should be similar to the simulation results. 
Thus effect (1) will also contribute little to the discrepancy 
on the final radial distribution of subhaloes. 

It is then reasonable to conclude that the over-predicted 
subhaloes at small radii is because either simulations still 
lack enough resolutions to resolve subhaloes in the central 
region of the host halo, or the model neglecting subhalo dis- 
ruption is not realistic. With respect to simulation, Springel 
et al. (2008) have shown that increasing the resolution does 
resolve more low-mass subhaloes, but the number of sub- 
haloes converges for given mass limit. Thus it is implausible 
that simulation resolution is responsible for this discrepancy. 
With respect to the model, defining a subhalo to be dis- 
rupted (or unbound) is very subjective, for example, most 
authors assume that subhalo is tidally disrupted if its mass 
is less than the initial mass within a radius f dial" a, but with 
a wide range of fdis between 0.1 ~ 2.0. As shown by Wetzel 
& White (2010), varying fdia has a huge impacts on the final 
radial distribution of subhaloes. 

In fact, there are more effects which can affect the abun- 
dance of subhaloes and their radial distribution, (i) Host 
halo formed in cosmological simulation always contains more 
than one subhalo, and the interaction between subhaloes will 
accelerate the disruption of subhaloes and reduce the num- 
ber of subhaloes at the inner host halo (e.g., Tormen et al. 
1998; Gnedin et al. 2004; Angulo et al. 2009). (u) Ludlow et 
al. (2009) have shown that small subhaloes are more likely 
to be ejected out to larger distances during the virialization 
of the host halo, thus producing a less concentrated distri- 
bution, (iii) Subhalo-subhalo mergers may be also effective 
to reduce the abundance of subhaloes (e.g., Kim et al. 2009), 
especially for the less massive subhaloes (e.g., Angulo et al. 
2009). Unfortunately, these processes are difficult to be in- 
cluded in the analytical model. 

Finally in this section we consider the dependence of 
subhalo radial distribution on their properties. Most N-body 
simulations have shown that the radial distribution of sub- 
haloes has no dependence on the present-mass of subhaloes 
(e.g., Gao et al. 2004b; Diemand et al. 2004; Springel et 
al. 2008), while others (e.g., De Lucia et al. 2004) found 
non-negligible dependence on subhalo mass. In Fig. [9] we 
show the fiducial model predictions with dependence on the 
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Figure 9. The dependence of subhalo radial distribution on subhalo properties. The left, middle and right panels show the dependence 
on the present-day mass, mass at accretion of subhaloes (both in unit of present host halo mass, as indicated in each panel) and their 
accretion redshift. In left panel, it also shows the distribution from the data of "Via Lactea" simulation as in grey lines. 



present-day mass (left panel), mass at accretion (middle 
panel) and accretion redshift (right panel). The subhaloes 
mass are in unit of host halo mass at z = 0, and are indi- 
cated in each panel. In the left panel, the grey lines are the 
results obtained from the public data of "Via Lactea" simu- 
latioijj, where only the dependence on subhalo present-day 
mass can be derived. 

The left panel shows that the radial distribution has a 
weak dependence on the present-day mass of subhalo. This 
is consistent with the results of De Lucia et al. (2004). In- 
terestingly, the "Via Lactea" simulation show an opposite 
trend that higher-mass subhaloes have a more concentrated 
distribution within r < 0.2Rvir- As "Via Lactea" simulated 
only one host halo, there is significant scatter on the radial 
distribution of subhaloes due to the limited number in given 
mass bins. When split by their initial mass (i.e., mass at 
accretion), the subhaloes have almost the same radial dis- 
tribution, as shown in the middle panel. This result also 
demonstrates that why the ridial distribution of suhaloes 
is independent of the dynamical processes (panel e and f 
of Fig. |8} , as the subhaloes catalogue is dominated by the 
less massive suhaloes. The right panel shows that there is a 
significant dependence on the age of subhaloes that the old 
population has a more concentrated distribution at small 
radius. This age- dependence was also found by Taylor & 
Babul (2005b), and they pointed out that this is because 
old subhaloes are accreted at lower distances when the host 
halo was smaller at high redshift. 



4 HALO OCCUPATION DISTRIBUTION OF 
SUBHALOES 

In order to further examine the validity of our model, 
we investigate the halo occupation distribution (HOD) of 
subhaloes. For this purpose, we produce subhalo popu- 
lation in a set of host haloes with mass M{z = 0) = 
10", 10^^ 10^^10", 10^'' ft- ^ Mo . For each host halo, we 
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Figure 10. The normalized second moment of the HODs of sub- 
haloes as a function of the host halo mass for three mass thresh- 
olds: Igm/M > -5, -4, -3. 



produce 100 realizations of the merger trees, with a resolu- 
tion of 10^ 10^ 10^ lO^", IO^Mq, respectively. The fiducial 
model M2 is used to model the evolution of each subhalo. 
We count the surviving subhaloes (A'^) in each realization 
and compute the quantity 



(N) 



1/2 



(16) 



http://www.ucolick.org/~diemand/vl 



for each host halo. The quantity a is the normalized sec- 
ond moment of subhalo's HOD (e.g., Kravtsov et al. 2004a; 
van den Bosch et al. 2005). For a Poissonian distribution, a 
should be unity, while distributions that are narrower (sub- 
Poissonian) or broader (super-Poissonian) have a < 1 and 
Q > 1, respectively. In the case of a « 1 (i.e., small devia- 
tion of a from unity), one should keep in mind that it does 
not necessarily indicate a Poissonian statistics (e.g., Boylan- 
Kolchin et al. 2010) and it needs further examination on the 
subhalo's HOD. Here we use the quantity a only for a consis- 
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tent comparison with the previous results (e.g., Z05; Zheng 
et aL 2005). 

In Fig. [To] we plot a as function of the mass of the 
host halo. We select subhaloes with different mass bins of 
Ig m/M > —5, —4, —3. Fig. 1 101 shows that a is close to unity 
independent of the mass threshold and the host halo mass. 
This result is inconsistent with the semi-analytical result of 
van den Bosch et al. (2005) and Z05, who found that the 
distributions of massive subhaloes in low-mass host haloes 
are significantly broader than Poissonian distribution. They 
pointed out that the discrepancy may be from the generic 
problem of conventional EPS formalism. However, our model 
predictions agree with the numerical results of Kravtsov et 
al. (2004a) and Zheng et al. (2005). We argue that this is 
due to the improvement of our model in two folds. Firstly, 
we employ a modified EPS formalism by Parkinson et al. 
(2008) which is shown to produce a well math to the merger 
history of haloes found from N-body simulations. Secondly, 
we assume that subhaloes never get disrupted but are can- 
nibalized with a time scale Tdf (Fig. [B]). It gives a better 
match to the observed distribution of satellite galaxies. 



5 DISCUSSION 

We have shown in Section 13.21 that the model predicts a 
more concentrated distribution of subhaloes than that seen 
from N-body simulations. This prediction is insensitive to 
the model assumptions for tidal stripping and dynamical pa- 
rameters. The same discrepancy was also obtained by Taylor 

6 Babul (2005), and they ascribed it to the numerical ef- 
fects of simulations. The recent high-resolution N-body sim- 
ulations (e.g., Springel et al. 2008) have shown that reso- 
lution is not the scapegoat for the low number density of 
subhaloes in the inner host halo. As pointed by Z05 that it 
is not the resolution but the model assumption that gives to 
the over-prediction. Z05 obtained a good match to the simu- 
lation result by implementing disruption of subhalo. Wetzel 
& White (2010) also have shown that decreasing the thresh- 
old of tidal disruption produces more subhaloes in the inner 
host halo. 

We believe that the good agreement between our model 
prediction and the observed distribution of the MW satel- 
lites is a coincidence. To get a better match to the distribu- 
tion of subhaloes from N-body simulations, we should firstly 
understand the importance of subhalo-subhalo interaction 
and subhalo-host interaction for the disruption of subhaloes, 
and more studies are needed to classify their contribution to 
this discrepancy. If these interactions are not enough to dis- 
solve subhaloes in the inner host halo, a more realistic model 
for subhalo disruption should be included in the model. Cur- 
rently, this is difficult to implement it. In one aspect, we need 
a more realistic and physical model for subhalo disruption, 
unlike those models to define disruption when the distance 
of subhalo to host center is less than some given radius be- 
cause it is arbitrary and dependent on the resolution of sim- 
ulations. On the other hand, subhalo disruption is easily to 
be confused with subhalo merging with central halo, which 
is often true for low-resolution simulations. For our exper- 
iment in this paper, we have to neglect the disruption and 
define subhalo merger using its angular momentum. This is 



adopted for making a fair comparison with the simulations 
of BK08. 

Though we have neglected subhalo disruption in our 
model, the result indicates that disruption is not important 
for real galaxies. This is seen from panel d of Fig|8] Indeed, 
Hydrodynamical simulations with baryon included have con- 
firmed that subhalo can survive the strong tidal disruption 
(e.g., Gnedin et al. 2004; Nagai & Kravtsov 2005; Weinberg 
et al. 2008; Dolag et al. 2009), and the predicted spatial dis- 
tribution of satellite galaxies is similar to that observed in 
galaxy clusters (e.g., Gao et al. 2004a). Actually this is not 
the only solution to this problem. Some have argued that ob- 
served satellites in the MW are biased tracers of subhaloes 
(e.g., Kravtsov et al. 2004b; Madau et al. 2008). The read- 
ers are referred to the review paper by Kravtsov (2010) for 
more discussions. 



6 CONCLUSION 

In this paper, we study the evolution of dark matter subhalo 
using an analytical model including simple descriptions for 
a few important processes, such as tidal stripping, dynam- 
ical friction, tidal heating. We tune the model parameters 
to fit the dynamical evolution of subhalo predicted by con- 
trolled N-body simulation. Then we combine these descrip- 
tions with merger trees from the EPS-based Monte-Carlo 
merger trees and study the subhalo population in a Milky- 
Way type halo, including the subhalo mass function and the 
radial distribution of subhaloes. 

Following Boylan-Kolchin et al. (2008), we define sub- 
halo to be merged with central halo when its angular mo- 
mentum reaches zero. We compare the predicted angular 
momentum evolution to the simulation results of BK08. We 
find that the mass loss of subhalo due to tidal stripping has 
great impact on its angular momentum evolution. A high 
tidal stripping efficiency. A, produces a fast decrease of sub- 
halo mass and a longer merger time scales. We further find 
that the mass of subhalo should not decrease continuously 
by tidal stripping, and better agreement with simulation can 
be obtained if the mass of subhalo keep fixed after two pas- 
sage of pericenter. We give a modified Coulomb logarithm 
using the fitting formula of BK08 to gain a well match to 
the subhalo merger time-scales. 

We compare the subhalo mass function to that from 
N-body simulations, and it is found that SHMF is mainly 
determined by the tidal stripping efficiency A, but depen- 
dent weakly on the Coulomb logarithm parameters. It is also 
insensitive to subhalo mass loss at its late stage of evolution, 
as the SHMF is dominated by recently accreted subhaloes. 
This is in good agreement with the results of van den Bosch 
et al. (2005) and Yang et al. (2009). 

The radial distribution of subhaloes is found to be more 
central concentrated than that from N-body simulations. 
This is a common prediction from analytical models if sub- 
halo disruption is not important (Taylor & Babul 2005b; 
Z05; Wetzel & White 2010). N-body simulations (e.g., Gao 
et al. 2004b; Springel et al. 2008) seems to indicate that 
resolution is not the scapegoat for the under-prediction of 
subhaloes at small radii. The interaction between subhalo 
and host halo is efficient to eject small subhaloes into the 
outer region of host halo (Ludlow et al. 2009) which can 
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partly resolves this discrepancy. To fully solve this problem, 
we need a more physical model for subhalo disruption and 
distinguish between disruption and merger, which are not 
feasible at the moment. 

We conclude that the better agreement between pre- 
dicted spatial distribution of subhaloes and observed satel- 
lites in the Milky- Way is a coincidence. But it also implies 
that real galaxy may not be tidally disrupted, as found from 
hydrodynamical simulations (e.g., Dolag et al. 2009) that 
condensation of baryon mass inside subhalo will increase the 
central density of subhalo, making it more resistant to tidal 
disruption. Another argument for the concentrated distribu- 
tion of satellite galaxies is that they are biased population of 
subhaloes found from simulations (e.g., Madau et al. 2008). 
A final solution to this issue should turn to the hydrody- 
namical simulation with high enough resolution and more 
realistic models for gas cooling and star formation in sub- 
haloes. 

Finally, we investigate the halo occupation distribution 
(HOD) of subhaloes with our improved model. The second 
moment of the HODs is close to unity, which disagrees with 
the results from semi-analytical model of van den Bosch et 
al. (2005) and Z05 but agrees with the HOD derived from 
numerical simulation by Kravtsov et al. (2004a) and Zheng 
et al. (2005). 
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